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ABSTRACT 



We present direct upper limits on continuous gravitational wave emission from 
the Vela pulsar using data from the Virgo detector's second science run. These 
upper limits have been obtained using three independent methods that assume 
the gravitational wave emission follows the radio timing. Two of the methods 
produce frequentist upper limits for an assumed known orientation of the star's 
spin axis and value of the wave polarization angle of, respectively, 1.9x10"^^ and 
2.2x10^^^, with 95% confidence. The third method, under the same hypothesis, 
produces a Bayesian upper limit of 2.1 x 10~^^, with 95% degree of belief. These 
limits are below the indirect spin-down limit of 3.3 x 10~^^ for the Vela pulsar, 
defined by the energy loss rate inferred from observed decrease in Vela's spin 
frequency, and correspond to a limit on the star ellipticity of ~ 10~^. Slightly less 
stringent results, but still well below the spin-down limit, are obtained assuming 
the star's spin axis inclination and the wave polarization angles are unknown. 

Subject headings: gravitational waves - pulsars: general 



Introduction 



We describe here a search for continuous gravitational radiation from the Vela pulsar 
(PSRB0833-45, PSR J0835-4510) in data from the Virgo detector VSR2 run, which began 
on 2009 July 7 and ended on 2010 January 8. Continuous gravitational waves (CW) can 
be emitted by a rotating neutron star through a variety of possible mechanisms, including 
non-axisymmetry of its mass distribution, giving rise to a time- varying quadrupole moment. 
Such emission would imply loss of rotational energy and decrease in spin frequency. Hence 
a pulsar's observed frequency spin-down can be used to place an indirect upper limit on its 
gravitational wave emission, named spin- down limit. While a recent search for CW radiation 
using LIGO data has been carried out for more than 100 known pulsars ( Abbott et al.||2010 ). 



the resulting upper limits have beaten the spin-down limit for only the Crab pulsar (Abbott 



et al. 2008 2010). A search over LIGO data for CW signals from the non-pulsing neutron 



star in the supernova remnant Cassiopeia A has established an upper limit on the signal 
amplitude over a wide range of frequencies which is below the indirect limit derived from 
energy conservation (Abadie et al. 2010). In this article we present upper limits on CW 



emission from the Vela pulsar that lie below its spin-down limit, making Vela only the 
second pulsar for which this experimental milestone has been achieved. The only previous 
targeted search for CW emission from the Vela pulsar was in CLIO data over the period 2007 
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February 12-28, which produced an upper hmit of ~ 5.3x10 several orders of magnitude 



above the spin-down hmit (Akutsu et ah 2008) 



Vela is observed to pulsate {frot — 11-19 Hz) in radio, optical. X-ray and 7-ray radiation 
and is associated with the Vela supernova remnant. The association of the pulsar to the 



supernova remnant was made in 1968 (Large et al. 1968 ) and was the first direct observational 
proof that supernovae can produce neutron stars. The Vela spin-down rate is frot — — 1.56x 
10~^^ Hz s~^, corresponding to a kinetic energy loss of Egd — 6.9x10^^ W, where the canonical 
value for a neutron star's moment of inertia, / = lO^^kgm^, has been assumed. This loss 
of energy is due to various mechanisms, including magnetic dipole radiation, acceleration 
of charged particles in the pulsar magnetosphere and possibly the emission of gravitational 
waves. In this analysis we assume a tri-axial neutron star rotating around a principal axis 
of inertia, so that the gravitational wave (GW) signal frequency is / = 2 frot (see Section [2]). 
With an estimated distance from the Earth of ~ 290 pc (Dodson et al. 2003), Vela is one of 



the nearest known pulsars. Assuming that all the observed spin-down is due to the emission 
of gravitational waves, we obtain the spin-down limit h^'^ = 3.29 x lO"^'^ for GW tensor 



amplitude at the Earth. With an estimated age of ~ llOOOyr (Caraveo & Bignami 1989) 



Vela is relatively young and could, in principle, have a significant residual non-axisymmetry 
from its formation. The spin-down limit on the signal amplitude can be converted into an 
upper limit on the star's equatorial ellipticity e (see Eq. 15 ). For Vela we have e'*'^ = 1.8x10"^. 



This value is far larger than the maximum allowed by standard equations of state for neutron 



star matter (Horowitz & Kadau 2009), but is comparable to the maximum value foreseen 



by some exotic equations of state (Owen 2005; Lin 2007 Haskell et al. 2007). Because 



of very effective seismic isolation (Acernese et al. 2010), Vela's GW emission frequency 
(/ ~ 22.38 Hz) is within the sensitive band of the Virgo detector; this frequency range is 
inaccessible to all other gravitational wave detectors to date. 

Vela is a particularly glitchy pulsar, with an average glitch rate of ~ 1/3 yr~^, making it 
important to know whether or not a glitch occurred during the VSR2 run. Vela is regularly 

Table 1. Position and estimated distance of Vela pulsar. For the position, parentheses 
give the la error on the final digit as produced by the TEMP02 fit; for the distance, the 
uncertainty estimated in (Dodson et al. 2003[ ) is quoted. Positional parameters refer to 

epoch (MJD) 54620. 



d [pc] 



08'' 35" 20.75438(3)" -45° 10' 32".9507(7) 287 (-17, +19) 
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monitored by both the Hobart radio telescope in Tasmania and the Hartebeesthoek radio 
telescope in South Africa. According to their observations, no glitch occurred during the time 
span of VSR2. Prior to VSR2 it last glitched on 2007 August 1, and it has since ghtched on 
2010 July 31 ( Buchner||2010 ). Observations from the Hobart and Hartebeesthoek telescopes 



have also been used to produce updated ephemerides for Vela, which are important given 
Vela's relatively large timing noise. If timing noise is a consequence of fluctuations in the 
star's rotation frequency, not taking it into account would result in an increasing mismatch 
over time between the signal and template phases, thus producing a sensitivity loss in a 
coherent search. In this search updated ephemerides have been computed using the pulsar 
software TEMP02 starting from the set of times of arrival (TOAs) of the electromagnetic 
pulses observed by the Hobart and Hartebeesthoek telescopes covering the whole duration 
of the VSR2 run. Including in the fitting process up to the second derivative of frequency is 
enough in order to have fiat post-fit residuals. The post-fit position and frequency parameters 
are shown in Table [T] and Table [2] respectively. The corresponding post-fit residuals rms 
amounts to a negligible 100 /is. 

Recent Chandra X-ray observations provide accurate determination of the orientation 
of the Vela spin axis. In ( Ng fc Romani||2008 ) estimates of the pulsar wind nebula's "position 



angle", ipp, and inclination bp are given: 

ijp = 130.63° ±0.05°, 

LP = 63.6° ±0.6°. (1) 

The "position angle" is related to the gravitational wave polarization angle (see Section 
[2]) by either ip = 180° + ipp or ip = ipp depending on the unknown spin direction. Our 
analyses are insensitive to rotations of ip by integer multiples of 90°, so the spin direction 
is not needed. The inclination angle calculated from the pulsar wind nebula Lp is taken to 
be the same as that of the pulsar l. The physics of pulsar wind nebulae is complex, and a 
model leading to the above fits has several uncertainties. Thus we perform separate searches 



Table 2. Spin frequency, spin-down rate and estimated age of Vela pulsar. Parentheses 
give the la error on the final digit of spin frequency and spin-down rate estimations as 
produced by the TEMP02 fit. Rotational parameters refer to epoch (MJD) 54620. The 
quoted precision is enough to determine the rotational phase to within about 0.012 cycles. 



frot [Hz] 


frot [HzS-1] 


frot [Hzs 2] 


Age [yr] 


11.19057302331(9) 


-1.5583876(4) X 10-" 


4.9069(9) X 10-22 


11000 
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for the GW signal from Vela, both assuming that the angles ip and l are known within the 
above uncertainties, and assuming that they are unknown. 

The remainder of this paper is organized as follows. In Section [2] we summarize the 
characteristics of the GW signals for which we search. In Section |3] we describe the data set 
used for the analysis. In Section |4] we briefly describe the three analysis methods used. In 
Section [5] we present the results of the analysis. In Section |6] we provide conclusions. Some 
more details on the analysis methods are given in the Appendices. 



2. The GW signal 

The continuous GW signal emitted by a triaxial neutron star rotating around a principal 
axis of inertia as seen from Earth is described by the following tensor metric perturbation: 

h(t) = /i+(t)e+ + /ix(t)ex, (2) 

where 

h+{t) = hoi jcos$(t) (3) 

hy^{t) = /lo cost sin <l>(t), (4) 



and e+ and ex are the two basis polarization tensors. They are defined, see e.g. (Misner 



et al. 1973), in terms of unit orthogonal vectors ex and ey where ex is along the x-axis of 
the wave frame, defined as the cross product s x n between the source spin direction s and 
the source direction n in the solar system barycenter (SSB). 

The angle l is the inclination of the star's rotation axis with respect to the line of sight 
and is the signal phase function, where t is the detector time, while the amplitude ho 
is given by 

ho = ^, (5) 

where I^z is the star moment of inertia with respect to the rotation axis, the equatorial 
ellipticity e is defined, in terms of principal moments of inertia, as e = h^^zJjDt^^ ^ jg 
star distance and / is the signal frequency. As the time-varying components of the mass 
quadrupole moment tensor are periodic with period half the star rotation period, it follows 
that / = 2frof 

The GW strain at the detector can be described as 

h{t) = h+{t)F+{t- ^) + h^ {t)F^ {t- ij), (6) 
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where the two beam-pattern functions, which are periodic functions of time with period of 
one sidereal day, are given by 

F+{t;^) = a{t) cos 2^ + b{t) sin 2tp (7) 
F^{t;^) = 6(t)cos2^-a(t)sin2^. (8) 

The two functions a(t),b(t) depend on the source position in the sky and on the detector 
position and orientation on the Earth. Their time dependency is sinusoidal and cosinusoidal 
with arguments t and 2 t, where f2® is the Earth angular rotation frequency; ip is the 
wave polarization angle defined as the angle from z x n to the x-axis of the wave frame, 
measured counterclockwise respect to fi, where z is the direction of the North celestial 



pole (see, e.g., the plot in (Prix & Krishnan 2009)). The effect of detector response on 



a monochromatic signal with angular frequency loq is to introduce an amplitude and phase 
modulation which determine a split of the signal power into five frequencies, uq, ijJo±Q^, loo± 
2fl^. The distribution of power among the five bands depends on the source and detector 
angular parameters. In Fig. [T] the power spectrum at the Virgo detector of an hypothetical 
monochromatic signal coming from the location of the Vela pulsar is shown for two assumed 
polarizations (pure "+" linear polarization and circular left handed polarization). 




22.38105 22.3811 22.38115 22.3812 22.38125 22.38105 22.3811 22.38115 22.3812 22.38125 

f[Hzl f[Hz] 



Fig. 1. — Power spectrum of an hypothetical monochromatic signal coming from the location 
of the Vela pulsar as seen from the Virgo detector. Left plot refers to a purely + signal; 
right plot to a circularly (left handed) polarized signal. 

To a very good approximation the SSB can be used as an inertial reference frame in 
which to define the signal phase. In this frame, with barycentric time T, the signal phase is 

<l>(T) = <l>o + 27r/o(T-To), (9) 

where the signal intrinsic frequency /o is a function of time due to the spin-down: 

/o(^) = /(°) + ^^(T-To)^ (10) 

n=l 
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where /'^"^ = ^j^\t=To- The time at the detector, t, differs from T due to the relative 
motion between the source and the detector and to some relativistic effects. Considering 



only isolated neutron stars we have the well-known relation (Lyne & Graham-Smith 1998 



Hobbs et al. 2006 Edwards et al. 2006) 



T = t + Ar + As + Ac 



:ii^ 



where 



A 



r ■ n 



R 



(12) 



is the classical Roemer delay, which gives the main contribution (r is the vector identifying 
the detector position in the SSB, while h is the unit vector toward the source). The term 
A^; is the Einstein delay which is the sum of two contributions, one due to the gravitational 
redshift produced by the Sun and the other due to the time dilation produced by the Earth's 
motion. As is the Shapiro delay due to the curvature of space-time near the Sun. Expressing 



the signal phase in the detector frame, by using Eq. 11, we can write the signal frequency at 
the detector as 



m 



1 d$(t) 
27r dt 



fo{t) 1 + 



V ■ n 



+ rel. corr. 



(l)r-n| 



(13) 



or smaller have been 



where v is the detector velocity vector and terms of order |/ 
omitted from the equation (though they are included in the analyses). 

A useful quantity to which to compare the upper limit on signal strength set in a given 



analysis is the so-called spin-down limit. It is computed (Abbott et al. 2007) assuming that 
all the observed spin-down is due to the emission of GW: 



/i^'^ = 8.06 X 10-19 Jssrfkpe 



|( Aot/Hzs-i) 

(/rot/Hz) 



(14) 



where Jag is the star's moment of inertia in units of 10^^ kgm^ and c/kpc is the star's distance 
from the Sun in kiloparsecs. It is an absolute upper limit to the amplitude of the GW signal 
that could be emitted by the star, where electromagnetic radiation is neglected. The spin- 
down limit on the signal amplitude corresponds to an upper limit on the star's ellipticity 
given by 



0.237 



10 



-24 



d 



kpc- 



(15) 



The Vela pulsar has a measured braking index n ~ 1.4 (Lyne et al 1996) and this, together 



with the estimation of its age, can be used to compute a stricter indirect hmit on the signal 



amplitude (Palomba 2000), which only holds under the assumption that the spin-down is 



due to the combination of emission of GW and magnetic dipole radiation, about 4 times 
lower than the spin-down limit. 
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Achieving sensitivity better than the spin-down hmit is an important milestone toward 
probing neutron star structure via gravitational waves. 



3. Instrumental performance in the VSR2 run 

We have analyzed calibrated strain data from the Virgo VSR2 run. This run (started 
in coincidence with the start of the LIGO S6 data run) began on 2009 July 7 21:00:00 UTC 
(GPS 931035615) and ended on 2010 January 8 22:00:01 UTC (GPS 947023216). The duty 
cycle was 80.4%, resulting in a total of ~ 149 days of science mode data, divided among 
379 segments. Science mode is a flag used to indicate when the interferometer is locked and 
freely running at its working point, with all the controls active and no human intervention. 
In Fig. |2] the fraction of total time covered by science data segments with duration longer 
than a given value is plotted. The longest segment lasts ~ 88 hours. 



1r 
0.8 

0.6 
0.4 
0.2 


10" 




10^ 10' 

duration [h] 



10^ 



Fig. 2. — Fraction of the total time covered by science data segments with duration larger 
than a given time. 



The detector showed a good sensitivity around the expected Vela signal frequency during 
the entire run. The sensitivity was typically within a factor of two of the target Virgo design 
sensitivity ( Accadia et a/.|[2010 ). Figure |3] shows the estimation of the power spectrum of the 
data, computed through an average of ~ 1, 000s periodograms after removal of some large 



outliers (see Section 4.3.1), on a 0.8 Hz frequency band around the expected frequency of 



the GW signal from the Vela pulsar for the entire VSR2 run. An instrumental disturbance 
right at the Vela signal frequency degraded the sensitivity by ~ 20% with respect to the 
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background. The source of this disturbance was seismic noise produced by the engine of the 
chiller pumps that circulate coolant fluid for the laser of the mirror thermal compensation 
system and it has been removed during the next Virgo VSR3 run. 




22.1 22.2 22.3 22.4 22.5 22.6 22.7 22.8 22.9 



f[Hz] 

Fig. 3. — Estimation of the power spectrum of VSR2 data in a 0.8Hz band around the 
expected Vela signal frequency. The expected signal frequency (vertical dashed line) is right 
in the middle of the frequency band affected by an instrumental disturbance, see text for 
more details. 



The data used in the analysis have been produced using the most up-to-date calibration 
parameters and reconstruction procedure. The associated systematic error amounts to 5.5% 

over the frequency range between 
~ 10 Hz and ~ 1 kHz, with lower uncertainties at the Vela frequency. The reconstructed 
data have a sampling rate of 20 000 Hz. However two more reconstructed data streams, 
sampled respectively at 16 384 Hz and 4096 Hz, were also produced to be consistent with 
LIGO/GEO sampling rates. 



in amplitude and ~ 50 mrad in phase (Accadia et a/.|2010 



The search methods 



Three different and largely independent analysis methods have been applied to this 
search: 1) a complex heterodyne method using Bayesian formalism and a Markov chain 
Monte Carlo (Abbott et al. 2010[ ); 2) a time-domain matched filter method using the J-'- 



statistic (Jaranowski et al. 1998) and a new extension known as the ^-statistic (Jaranowski 



& Krolak 2010); and 3) a matched filter method applied to the signal's Fourier components 
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at the five frequencies to wliicli tlie signal is spread by tlie sidereal modulation (Astone et al. 



2010). 



There are several reasons to use different methods in the search for CW signals, provided 
they have comparable performance. First, it makes it easier to cross-check each method by 
comparing the analysis outputs, even at intermediate steps. Second, different methods can 
be more suitable, or efficient, for given characteristics of the data to be analyzed, or for 
given characteristics of the signal emitted by a source; e.g. a method can be more robust 
against noise non-stationarity with respect to another. Third, in case of detection with a 
given analysis it will be of paramount importance to confirm the detection with one or more 
independent analyses. 

In the analyses described in this paper we observe consistent results from the three 
methods, which provide valuable cross-checks. 

All the analyses clean the data in some way to remove large transient outliers. This is 
necessary, as large short-duration transients will skew noise estimates and adversely affect 
results. The amount of data removed during cleaning is negligible compared to the total 
data span and would produce a decrease of the signal-to-noise ratio of a signal present in 
the data of less than 1%. 

Among the three methods two different approaches have been used towards setting up- 
per limits. In the heterodyne method the posterior probability for the signal parameters is 
calculated, from which degree of belief (or credibility) regions can be set to give limits on 
particular parameters {e.g., an upper limit on ho can be set by finding the value that bounds 
a given percentage of the probability). In the two other analyses a frequentist approach is 
used and upper limits are set through Monte Carlo methods where many simulated signals 
with different amplitude and randomly varying parameters and frequency near the expected 
one from the Vela are added to the data. These two approaches should produce quantita- 



tively similar results, see for example (Abbott et al. 2004), but they are answering different 



questions and therefore cannot be meaningfully combined. 

The three analysis methods are described in the following sections of this paper. 



4.1. Complex heterodyne 



This method, developed in (Dupuis &; Woan 2005), provides a way to reduce the search 
dataset to a manageable size, and use it to perform Bayesian parameter estimation of the 
unknown signal parameters. 
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4- 1.1. Data reduction 



The known signal phase evolution (Eq. |9]) is used to heterodyne the data, changing the 
time series detector data x{t) = h{t) + n(t), where h{t) is the signal given by Eq. [6] and n{t) 
is the noise, to 

x'it) = x(t)e-'[*W-*ol, (16) 

giving a complex dataset in which the signal is given by h'(t) = /i(t)e~*t**^*^~*°l. The now 
complex signal is 

h'(t) = /lo ^^-FV(1 + cos^i) cos$o + ^-P'x cosisin$o^ + 

iho ^-F+(l + cos^i) sin$o — -F^ cost cos , (17) 

where and Fx are given by Eqs. [7] and [8j This heterodyne therefore removes the fast- 
varying part of the signal (the time dependent part of Eq. [9]) leaving a complex data stream 
with the signal shifted to zero frequency (setting aside small offsets due to the diurnal ampli- 
tude modulation of the signal from the detector beam pattern). In practice this heterodyne 
is performed in a two stage process. First a coarse heterodyne is performed using the phase 
evolution calculated assuming a stationary frame. This data is then low-pass filtered (in 
this case using a 9th order Butterworth filter with a 0.25 Hz knee frequency) and heavily 
downsampled from the original rate of 16 384 Hz to 1 Hz. A second stage of heterodyne takes 
into account the signal's modulation due to the Earth's motion and relativistic effects (see 



Eq. 13). The data are then further downsampled from 1 Hz to 1/60 Hz by taking the mean 



of 60 samples, which has the effect of an additional low pass filter. 



4-1.2. Data cleaning 

The fully heterodyned data are cleaned to remove the largest outliers, by discarding 
points with absolute values greater than 5 times the standard deviation of the data. This 
cleaning is performed twice to combat the effect of extreme outliers (many order of magnitude 
larger than normal) skewing the standard deviation estimate. This removes ~ 0.05% of the 
data. 



For the parameter estimation, as in dAbbott et al.|2007|[20T0| , the hkelihood calculation 



assumes the data is stationary for contiguous 30 minute segments, although shorter segments 
of 5 minutes or more are also included to account for shorter stretches of data at the end of 
longer contiguous segments. This contiguity requirement removes a further ~ 0.2% of the 
heterodyned data, which is within segments shorter than 5 minutes. 
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4.1.3. Parameter estimation and upper limits 

This new, and far smaller, 1/60 Hz-sampled dataset is then used to estimate the four 
unknown signal parameters ho, $0) cos 6, and ip. These are estimated using a Bayesian 
formalism, with a Students-i-like distribution for the likelihood (formed by marginalizing a 
Gaussian likelihood over an unknown noise standard deviation) given the heterodyned data 



and a signal model from Eq. 17, and specific priors (see below) on these parameters. This 



posterior probability volume is explored using a Markov chain Monte Carlo (Abbott et al. 



2010), which gives posterior probability distribution functions (PDFs) on each parameter 



marginalized over the three others. 

In this analysis two different sets of independent priors are used for the parameters. 
In one case uniform priors on all four parameters are set - for the angular parameters this 
means that they are uniform across their allowable ranges, but for the lower bound is 
zero, and the upper bound is set at a level well above any values that could be consistent 
with the data. For reasons set out in Section [l] the other case sets the priors on ip and cos l 
to be Gaussians given by Eq. [ij whilst keeping the ho and $o priors as uniform. 

The marginalised ho posterior, p{ho\d,I), can be used to set an upper limit in the 
amplitude by finding the value of /iq' that bounds (from zero) the cumulative probability to 
a given degree-of-belief, B, 

B= / p{ho\d,I)dho. (18) 
Jo 

Here we set 95% degree-of-belief upper limits. Due to the fact that the MCMC is finite in 
length there will be small statistical uncertainties between different MCMC runs, which for 
cleaned data we find to be < Ix 10^^^. The difference in results between using cleaned and 
non-cleaned data, as above, is within the statistical uncertainty from the MCMC. 



4.2. and Q statistics method 



The second search method uses the J-' and Q statistics developed in (Jaranowski et al. 



1998) and (Jaranowski & Krolak 2010). These statistics are used to perform maximum- 



likelihood estimation of signal parameters and to obtain frequentist upper limits on the 
signal amplitude. 
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4-2.1. Data reduction 



The description of how to compute the J-" and Q statistics from time-domain data is 
given in ( Jaranowski et aL]|1998 ) and ( Jaranowski &: Kr61ak|2010 ). The J-' statistic is apphed 
when the four parameters ho, $o; and i are assumed to be unknown. When the orientation 
of the spin axis of the Vela pulsar and the wave polarization angle are known and given by 
Eqs. [T} the Q statistic is used instead. 

We have refined the application of these statistics to account for two features of the 
current search. Firstly, the VSR2 data that we analyze are not stationary (see Figure |6]), so 
the statistics must be adjusted to de-emphasize noisy periods. Secondly, we use as our input 



data the complex-valued coarse heterodyne data described in Section |4.1[ so the statistics 
must be generalized to deal with complex data. These effects can be taken into account in 
J-" and Q statistics formalism in a straightforward way derived explicitly in Appendix O 



resulting in the generalized forms of the J-' and Q statistics given by Eqs. Cll and C15 



respectively. These generalized forms of the statistics are used to search VSR2 data for a 
gravitational wave signal from the Vela pulsar. 



4.2.2. Data cleaning 



The coarse heterodyne data that we analyze with the J-" and Q statistics contains a 
small number of outliers that must be discarded. To identify these outliers we have used an 



iterative method called the Grubbs' test (Grubbs 1969) explained in detail in Appendix D 



Application of the Grubbs' test resulted in removal of 0.1% of the total data points in input 
data, amounting to a negligible loss of signal-to-noise ratio of any continuous signal present 
in the data. 



4.2.3. Parameter estimation and upper limits 

In the frequentist approach, a signal is detected in the data if the value of the J-" or ^ 
statistic exceeds some threshold corresponding to an accepted false alarm probability (1% 
in this analysis). When the values of the statistics are not statistically significant, we can 
set upper limits on the amplitude of the GW signal. We choose a frequentist framework 
by computing the amplitude of a signal that, if truly present in the data, would produce 
a value of the detection statistic that in 95% of the cases would be larger than the value 
actually found in the analysis. To obtain the upper limits on Hq we follow a Monte Carlo 



method described in (Abbott et al. 2004). That is, we add simulated GW signals to the VSR2 



- 20 - 



data and determine the resulting values of the statistics. The parameters of the simulated 
signals are exactly the same as for Vela, except for the gravitational-wave frequency which is 
randomly offset from twice the Vela spin frequency. For the T statistic case, the parameters 
i\} and cost are chosen from a uniform distribution, whereas for Q statistic case they are 
fixed to the values estimated from X-ray observations (see Eqjl]). We calculate the upper 
limits corresponding to the obtained values of the statistics by interpolating results of the 
simulation to find the value for which 95% of the signals have a louder T- or ^-statistic 
value than that obtained in the search. To estimate the statistical errors in the upper limits 
from the Monte Carlo simulations we have followed the method presented in Section IVE 



of (Abbott et al. 2004) by performing an additional set of injections for the amplitude 



around the obtained upper limits. 

In the case that a statistically significant signal is detected we can estimate unknown 
signal parameters. In the case of the T statistic search the maximum likelihood estimators 
of the amplitudes are obtained by equations |C12[ These amplitude estimates are then 
transformed into estimates of parameters /iq, $05 ^'iid i using Eqs. 23 of (Jaranowski & 



Kr61ak|[20T0| . In the case of the Q statistic search, where parameters i^) and i are assumed to 



be known, the amplitude estimator is obtained by Eq. C17 and estimates of the parameters 



/lo and $0 are calculated from Eqs. 7 of (Jaranowski & Krolak 2010) 



4.3. Matched filter on the signal Fourier components 

The third search method uses the Fourier amplitudes computed at five frequencies where 
the signal would appear due to sidereal amplitude modulation, and applies a matched filter to 
this 5-point complex data vector. Further details can be found in Appendix |A] and in (Astone 
et al.||2010l). 



^.5". i. Data reduction 

The starting point for this method is a short Fourier transform database (SFDB) built 
from calibrated strain data sampled at 4096 Hz ( Astone et al.|[2005 ). The FFTs have a dura- 
tion of 1024 s and are interlaced by 50% and windowed with a flat top - cosine edges window. 
From the SFDB a small band (0.2 Hz in this analysis) around the frequency of interest is 
extracted from each FFT. The SFDB contains, among other information, the position and 
the velocity of the detector in the SSB at the center time of each FFT. Each frequency do- 
main chunk is zero-padded and inversely Fourier-transformed to obtain a complex time series 
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with the same samphng time of the original time series, but with a spectrum different from 



zero only in the selected band {i.e., it is an analytical signal, see e.g., (Astone et al. 2002)). 
Then, for each sample, the detector position in the SSB is computed, by interpolating with 
a S""*^ degree polynomial. The Doppler and Einstein effects can be seen as a varying time 
delay A(t). A new non- uniformly-sampled time variable t' with samples t'^ = ti + A{ti) is 
computed. The spin-down is corrected by multiplying each data chunk by e"*^'^"'* where 
^(f>sdit') = 27r + Then the data are resampled at equal intervals in t'. The 

final complex time series has a sampling frequency of IHz. At this point, a true GW sig- 
nal would be sinusoidal with a sidereally-modulated amplitude and phase, as described in 
Sec.Q, containing power at the nominal source frequency and in lower and upper sidebands 
of ±^0, ±2^0. The Fourier coefficients at these five frequencies are taken to form a complex 
data 5- vector X. 

The detection method described here relies on a description of the GW signal given in 
Appendix |X| When the polarization angle ip and the inclination angle of the star rotation 
axis L are unknown, we use a procedure that we denote four-degrees of freedom detection, 
in which the two signal 5- vectors A'^,A^, corresponding to the + and x polarizations and 
defined in Appendix |A| are numerically computed and projected onto the data 5- vector X: 

X-A+ 

= -JJ+f (19) 
X ■ 

H. = -jj^- (20) 

The output of the two matched filters are the estimators of the amplitudes Hqc^'^^Hj^, Hqc^'^^H-^ 
The final detection statistic is defined by 

S= |A+n^+|' + |A><Hi^x|'. (21) 



More details can be found in (Astone et al. 2010). 



If estimations of i}) and l provided by X-ray observations (Section [T]) are used, we can 
apply a simpler procedure that we call a two-degrees of freedom detection. In this case the 
signal is completely known, apart from an overall complex amplitude H = Hqc^'^^. Then, 
the template consists of just one 5- vector A = H^A^ + Hy.A^, where H^, H^. are given 
by Eqs. |A3[ and only one matched filter must be applied to the data 5- vector X: 

^-W' (22) 

which provides an estimation of the signal complex amplitude. The detection statistic is 
then given by S = \H\'^. 
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4-3.2. Data cleaning 

In addition, various cleaning steps were applied to the data. The data can be modeled 
as a Gaussian process, with slowly varying variance, plus some unmodeled pulses affecting 
the tails of data distribution. The cleaning procedure consists of two parts. First, before the 
construction of the SFDB, high-frequency time domain events are identified after applying 
to the data a first-order Butterworth high-pass bilateral filter, with a cut-off frequency of 
100 Hz. These events are then subtracted from the original time series. In this way we do 
not reduce the observation time because we are simply removing from the data the high- 
frequency noisy component. The effect of this kind of cleaning has been studied in data from 
Virgo Commissioning and Weekly Science runs and typically reduces the overall noise level 
by up to 10 — 15%, depending on the quality of the data ( Acernese et a/.|2009 ). After Doppler 



and spin-down correction, further outliers that appear in the small band to be analyzed are 
also removed from the dataset by using a threshold of ±5xl0~^^ on the data strain amplitude, 
reducing the amount of data by ~ 1.3%. Slow non-stationarity of the noise is taken into 
account by applying a Wiener filter to the data, in which we estimate the variance of the 
Gaussian process over periods of ~ 1 000 s, and weight the data with its inverse in order to 
de-emphasize the more disturbed periods. 



4.3.3. Parameter estimation and upper limits 

Following the frequentist prescription, the value of S obtained from the search is com- 
pared with a threshold S* corresponding to a given false alarm probability (1% in this 
analysis). If S > S*, then one has a potential signal detection deserving deeper study. In 
the case of signal detection, the signal parameters can be estimated from H^, H^., using the 
relations shown in Appendix [B} If the measured S value lies below the threshold, we can set 
an upper limit on the amplitude of a possible signal present. 

The determination of upper limits is carried out via Monte Carlo simulations similar 



to the limit determination described in Section 4.2.3 In the case of 4 degrees of freedom 
(d.o.f.), the unknown parameters, ip and cost were taken to be uniformly distributed. The 
analysis method allows us to establish an upper limit for the wave amplitude Hq defined 
in Appendix |A| This was translated into an upper limit on ho, under the assumption that 



the source is a tri-axial neutron star, using Eq. A5 after maximising the factor under the 
square root with respect to the inclination angle. In this way the upper limit we obtain is 
conservative. In the 2 d.o.f case we compute the upper limit by using for ip and l the values 
given in Eqs. [T] 
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The statistical error associated with the Monte Carlo simulations is estimated as half of 
the difference between the two signal amplitudes that bound the 95% confidence level. The 
grid in the amplitude of the injected signals has been chosen fine enough that the resulting 
statistical error is about one order of magnitude smaller than the systematic error coming 
from calibration and actuation uncertainty. 



5. Results from the searches 

In the analyses all available science mode data recorded by Virgo were used. No ev- 
idence for a continuous gravitational wave signal was seen using any of the three analysis 
methods described in Section |4| We have therefore used the data to set upper hmits on the 
gravitational wave amplitude. 

For the complex heterodyne method (Section |4.1[ ) the marginalised posteriors for the 
four parameters, using the two different priors, are shown in Figs. |4] and [5j The presence 
of a detectable signal would show up as a posterior distribution in ho that is peaked away 
from ho = 0. The observed distributions are consistent with no signal being present. The 
95% credible limits on ho are shown and have values 2.4 x 10~^^ and 2.1 x 10~^^ respectively 
(note that the strongly-peaked distributions for cost and ijj in Fig. |4] are simply the restricted 
priors placed on those parameters). 



For the J-" and Q statistics (Section 4.2), the values obtained were consistent with false 
alarm probabihties of 22% and 35%, respectively. Since these probabilities are far above our 
1% false alarm threshold, we conclude that the data are consistent with the absence of a 



signal. Using the Monte Carlo method described in Sec. |4.2.3| we set 95% confidence upper 
limits on ho of 2.4x10"^^ and 2.2x10"^'^, respectively. 



For the matched filter on Fourier components (Section 4.3), the values computed for the 
4 d.o.f and 2 do./ statistics were consistent with false alarm probabilitis of 46% and 40%, 
respectively. Again we conclude that the data are consistent with the absence of a signal. 
We obtain 95% confidence upper limits on ho of 2.2 x 10~^^ and 1.9 x 10~^^, respectively. 

The results for all three analyses are summarized in Table |3} which also includes the 
systematic uncertainty in the upper limit from calibration and actuation uncertainties. For 
each analysis, results are given both for the case in which ip and cos l are assumed to be 
known [i.e. with restricted priors) and unknown {i.e. with unrestricted priors). 



We emphasize once again that the two results for the complex heterodyne method are 
Bayesian 95% credible hmits on ho, while the Q, J-", 2 d.o.f, and 4 d.o.f results are frequentist 
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95% confidence upper limits. While we would expect the two types of upper limit to be 

similar in value, they are not directly comparable, because they address different questions. 
The Bayesian question asks: "Given our priors and our data, for what value of ho are we 95% 
certain that any true signal lies below that value?" The frequentist question asks: "Above 
what value of ho would a signal produce a larger value of our statistic 95% of the time?" 
The subtle difference between these questions means that they may give different answers 
for the same data, and we should not read too much into the fact that in this search the two 
approaches gave very similar numbers. 
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Fig. 4. — The posterior PDFs for the pulsar parameters ho, $oi cost and ijj for 
PSR J0835— 4510, produced using restricted priors on cost and ip with the complex het- 
erodyne method. The vertical dashed line shows the 95% upper limit on ho. 
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Fig. 5. — The posterior PDFs for the pulsar parameters ho, cost and -0 for 

PSR J0835— 4510, produced using uniform priors for cost and ijj across the range of their 
possible values with the complex heterodyne method. The vertical dashed line shows the 
95% upper limit on ho. 

5.1. Validation with hardware injections 

All three pipelines used in the analysis have been tested with both software and hardware 

injections of CW signals in the VSR2 data. In particular wc discuss here hardware injections. 
For the entire duration of the run 13 CW signals (named PulsarO-12) have been injected 
in the Virgo detector by sending the appropriate excitations to the coils used to control 
one mirror's position. These signals were characterized by various amplitudes, spanned a 
frequency range from ~ 20 Hz to ~ 1400 Hz, and covered a range of values for the spin-down 
/ from ~ —4x10"-'^^ Hzs"-^ to ~ —2.5x10"^ Hzs~^. The corresponding source position (a, 5), 
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inclination t of the source spin axis, and polarization angle tp were chosen randomly. All 
the injected signals have been generated using the same software as the signals injected in 
LIGO S5 and previous runs. Injected signals PulsarO-9 have also the same parameters as the 
LIGO injections, while PulsarlO-12 have very low frequency and have been injected in Virgo 
only. The three pipelines were exercised on several of these simulated signals. The pipelines 
have been able to detect the signals and to estimate their parameters with good accuracy 
when the signal-to-noise ratio (SNR) is sufficient. In particular, in Tables |5|j7] we report 
the results obtained for PulsarS, characterized by a very small spin-down and high SNR, 
PulsarS with low frequency, very small spin-down and relatively low SNR and PulsarS with 
high spin-down and SNR. The frequency parameters for these three injections are given in 
Table |4j There is good agreement between the true and recovered signal parameters. With 
the method based on matched filtering on the signal Fourier components the estimation of 
the signal absolute phase is not straightforward. 



Conclusions 



In this paper we present the results of the analysis of Virgo VSR2 run data for the search 
of continuous GW signals from the Vela pulsar. The data have been analyzed using three 
largely independent methods and assuming that the gravitational wave emission follows the 
radio timing. For an assumed known orientation of the star's spin axis and value of the 
polarization angle, two methods have determined frequentist upper limits at 95% confidence 
level of, respectively, 1.9x10"^^ and 2.2x10"^^. The third method has determined a Bayesian 
95% degree-of-belief upper limit of 2.1 x 10"^'^. The lowest of these is about 41% below the 
indirect spin-down limit. It corresponds to a limit on the star ellipticity of 1.1 x 10^'^, 
which is well above the maximum equatorial ellipticity that a neutron star with a 'standard' 
equation of state can sustain, but comparable to the maximum value permitted by some 



exotic equations of state (Owen 2005; Lin 2007 Haskell et al. 2007). Given that the power 
emitted in GW is Egw = -^^I^e'^f, our results constrain the fraction of spin-down 
energy due to the emission of GW to be below 35%. For an unknown orientation of the star's 
spin axis and polarization angle the two frequentist upper limits are, respectively, 2.2x10"^^ 
and 2.4x10"^^ while the Bayesian upper limit is 2.4x10"^^. The lowest of these is about 33% 
below the spin-down limit. In this case the limit on the star ellipticity is 1.2 x 10"'^, while 
the corresponding limit on the fraction of spin-down energy emitted through GW is 45%. 
These numbers assume the canonical value for the star moment of inertia, / = lO^'^kgm^. 
However, the theoretically predicted values of I vary in the range ~ 1 — 3 x 10'^'^ kg m^ 



(Abbott et al. 2010), so our upper limit on the ellipticity can be considered as conservative. 
Such ellipticities could also be sustained by internal toroidal magnetic fields of order 10^^ G, 
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Table 3. Estimated 95% upper limit on Hq for PSR J0835— 4510 from the three different 
analysis methods (the horizontal line separates Bayesian from frequentist results). The 
systematic error on amplitude from calibration and actuation amounts to ~ 5.5%, as 
discussed in Section |3| This corresponds to an uncertainty on the upper limits of about 
±0.1 X 10~^^. For all upper limits the statistical error, associated with the Monte Carlo 
simulations used to the establish the limit itself, is about one order of magnitude smaller. 



Analysis method 


95% upper limit for 


Heterodyne, restricted priors 


(2.1 ±0.1) X 10- 


-24 


Heterodyne, unrestricted priors 


(2.4 ± 0.1) X 10' 


-24 


^-statistic 


(2.2 ±0.1) X 10- 


-24 


J^-statistic 


(2.4 ±0.1) X 10- 


-24 


MF on signal Fourier components, 2 d.o.f. 


(1.9 ±0.1) X 10- 


-24 


MF on signal Fourier components, 4 d.o.f. 


(2.2 ±0.1) X 10- 


-24 



Table 4. Frequency and positional parameters for the hardware injections (/ = for all 
the injections). The reference time epoch for the source frequency is MJD=52944 for all 
the injections. The optimal signal-to-noise ratio (SNR) is also given. 



Name 


f [Hz] 


/ [Hzs-i] 




a [deg] 


6 [deg] 


SNR 


PulsarS 


108.8571594 


-1.46x10" 


-17 


178.372574 


-33.436602 


192 


Pulsarb 


52.80832436 


-4.03x10- 


18 


302.626641 


-83.8391399 


40 


PulsarS 


194.3083185 


-8.65x10- 


-9 


351.389582 


-33.4185168 


197 



Table 5. Estimated parameters for hardware injection PulsarS from the three different 

analysis methods. 



Method 


ho. found 
ho.inj 


t [l^inj = 1.651] 


V' [V'mj = 0.444] 


^-0 ["^0,1,13 = 5.53] 


Heterodyne 


0.97 


1.67 


0.43 


5.55 


J^-statistic 


0.96 


1.65 


0.44 


5.54 


MF on signal Fourier comp., 4 d.o.f. 


0.96 


1.66 


0.44 


* 
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depending on the field configuration, equation of state, and superconductivity of tfie star 



(Akgun fc Wassermann 2007 Haskell et al. 2008 Colaiuda et al. 2008 Ciolfi et al. 2010). 



Then, our results have constrained the internal toroidal magnetic field of the Vela to be less 
than of the order of that value (it must be stressed, however, that the stability of a star with 
an internal field much larger than the external one is still an open issue). Vela is the second 
young pulsar for which the spin-down limit has now been beaten. 

A more stringent constraint on the emission of GW from the Vela pulsar may be estab- 
lished by analyzing data of the next Virgo+ run (VSR4) which is tentatively scheduled for 
summer 2011 and should last a few months. This run, assuming the planned sensitivity is 
reached, could be able to probe values of the Vela pulsar ellipticity below a few units in 10~^, 
corresponding to a fraction of spin-down energy emitted through the emission of GW below 
a few percent. We note that this run will also provide interesting results for several other 
low frequency pulsars. In particular, it could allow detection of GW from the Crab pulsar 
and J1952+3252 if their ellipticities are larger than ~ 10~^, a value nearly compatible with 
the maximum deformation allowed by standard neutron star equations of state. 

Second-generation detectors are expected to have a still better sensitivity at low fre- 



quency. Advanced Virgo (Acernese et al. 2009) and Advanced LIGO (Harry & the LIGO 



Scientific Collaboration 2010), which should enter into operation around 2014-2015, in one 
year could detect a GW signal from the Vela pulsar if its ellipticity is larger than a few times 
10~^, the corresponding fraction of spin-down energy emitted through GW being below a 
few times 10~^ in this case. 

The possibility of building a third generation GW detector, with a sensitivity a factor of 
10 or more better than Advanced detectors in a wide frequency range, is also being studied. 
The Einstein Telescope ( Punturo et al.||2010 ), which is currently at the stage of design study, 
is expected to release its first science data around 2025-2027. It should be able to detect 
GWs from the Vela pulsar, using one year of data, for ellipticity larger than 4x10^^ — 10~^, 
depending on the detector configuration that will be chosen. 



Table 6. Estimated parameters for hardware injection PulsarS from the three different 

analysis methods. 



Method 


found 


i iv) Wnj = 1.089] 


^ [Anj = -0.364] 


$0 [*0,z,y = 2.23] 


Heterodyne 


0.90 


0.99 


-0.27 


2.05 


J^-statistic 


0.89 


0.98 


-0.27 


2.10 


MF on signal Fourier comp., 4 d.o.f. 


0.97 


0.96 


-0.26 


* 
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A. An alternative formalism to describe a continuous GW signal 

The continuous GW signal emitted by a generic rotating rigid star can be described 
by a polarization ellipse. The polarization ellipse is characterized by the ratio = | of its 
semi-minor to its semi-major axis and by the angle i/j defining the direction of the major 
axis. The angle tjj is the same introduced in Sec. [2] The ratio r] varies in the range [—1, 1], 



Table 7. Estimated parameters for hardware injection PulsarS from the three different 

analysis methods. 



Method 


ha. found 


i (ry) [t,„, = 1.497] 


^ [Anj - 0.170] 


$0 [*0,,:nj = 5.89] 


Heterodyne 


0.97 


1.49 


0.18 


5.90 


J^-statistic 


0.95 


1.49 


0.17 


6.07 


MF on signal Fourier comp., 4 d.o.f. 


0.98 


1.50 


0.17 
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where 77 = for a linearly polarized wave and 17 = ±1 for a circularly polarized wave. The 
(complex) signal can be expressed as 

h(t) =/7o(^+e+ + i/xex)e^*(*), (Al) 

where e_|_ and ex are the two polarization tensors and the plus and cross amplitudes are 
given by 

cos 2iIj -ir] sin 2iIj 
+ " ^ ' 

sin 2tlj + ir] cos 2ij 
= /— ■ (A3) 

If we consider, as in Section [2| a triaxial neutron star rotating around a principal axis of 
inertia the following relations among Hq, t] and ho, t hold: 

2 cos i , ^ ^ , 

1 + cos 

Ho = yVl + Gcos^i + cos^i. (A5) 
In terms of + and x components we have 

^+ 7 (l + cos^i) 
^{Hy,^^=o) = -jr = FT — • (A7) 

-no -no 

In this formalism the complex gravitational strain at the detector is given by 

h{t) = Ho (A+(t)/7+ + Ax(t)i/x) e**^*\ (A8) 

where 

= F+(V^ = 0) (A9) 
Ax=Fx(V' = 0). (AlO) 



After Doppler and spin-down corrections, as described in Sec. |4.3[ we have: 

h{t) = Ho iA4t)H+ + Ax(t)ifx) e*("°*+*°^ (All) 

We now introduce the signal 5-vectors for the + and x components, A+, , given by the 
Fourier components, at the 5 frequencies produced by the amplitude and phase modulation, 
of the detector response functions A^, A^. It is straightforward to see that the signal in the 
antenna is completely defined by the 5-components complex vector 

A = i7oe^*° {H+A+ + H^A"") . (A12) 



More details can be found in (Astone et al. 2010). 
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B. Pcirameter estimators for MF on signal Fourier components 

Once the two estimators have been computed from the data, if a detection is 

claimed, the signal parameters Hq, t], ip can be estimated using the following relations. The 
estimator of the signal amplitude is given by 



= + (Bl) 



Introducing the quantities 



H+-Hy,^A + iB, (B2) 
l^+l^-l^xT^C, (B3) 

where the scalar product is between two complex numbers and includes a complex conjuga- 
tion of one, the estimation of the ratio between the axes of the polarization ellipse is 

. -1 + VI - 

V = ^ , (B4) 

while the estimation of the polarization angle can be obtained from 

cos (4^^) = , ^ (B5) 
2A 

sin (4^;) = , (B6) 



C. J-" and Q statistics for complex heterodyne data in non-stationary, 

uncorrelated noise 

We assume that the noise in the data is Gaussian and uncorrelated. In order to take 

into account non-stationarity of the data, wc assTime that each noise sample n{l) in the data 
time series is drawn from a Gaussian distribution with a variance ct^(/). We assume that 
the Gaussian distributions in question have zero means. Thus the autocorrelation function 
K{1, 1') for the noise is given by 

K{l,l') = a\l)6a', (CI) 

where /, /' are integers and where 6ui is Kronecker's delta function. Let us first assume that 
the signal h{l) is completely known and that the noise is additive. Thus when the signal is 
present the data take the following form 



x{l)^n{l) + h{l). 



(C2) 
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For Gaussian noise the optimal filter q{l) is the solution of the following (integral) equation 
(see (Jaranowski & Krolak 2009), p. 72) 

N 

h{l) = Y,K{lX)q{l'), (C3) 
v=i 

where is the number of data points. Consequently we have the following equation for the 
filter q{l) 

^(0 = ^ (C4) 
a\l) 

and the following expression for the log likelihood ratio In A 

\nK[x\ = {hx) -]^{h^), (C5) 



where the operator (■) is defined as 



Thus we see that for non-stationary Gaussian noise with the autocorrelation function (CI) 
the optimal processing is identical to matched filtering for a known signal in stationary 
Gaussian noise, except that we divide both the data and the filter by time-varying standard 
deviation of the noise. This may be thought as a special case of whitening the data and then 
correlating it using a whitened filter. The method is essentially the same as the Wiener filter 
introduced in Section 4^ The generalization to the case of signal with unknown parameters 
is immediate. 

In the analysis we use complex heterodyne data Xhet^ 

a:,et(0=a:(/)e-^*-('), (C7) 

where (^het is the heterodyne phase {^het can be an arbitrary real function). Thus we rewrite 
the J-' and Q statistics and amphtude parameter estimators using complex quantities. We 
introduce complex amplitudes Aa and A^, 

A, = A + zAs, (C8) 
A = A2 + iAi, (C9) 

where the amplitudes A^, = 1, 2, 3, 4 are defined by Eqs. 23 of ( Jaranowski fc Kr61ak|[20l0 ) 
and we also introduce the complex filters 

K{1) = a(/)e-'[*«-*'-*™ (CIO) 
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where a and b are amplitude modulation functions (see Eqs.[7]and[8]) defined by Eqs. 12 and 
13 in (Jaranowski et al. 1998), and $(/) is the phase defined by Eq. M 



(Cll) 



The J-'-statistic takes the following form 

J, ^ {b'^)\{Xhet ha)\^ + {a^)\{Xhet h)\'^ - 2{ab)^{{Xhet ha){Xhet h)*] 

{a^){b^) - {abf 

and the complex amplitude parameter estimators are given by 

^ (6^) {Xhet ha)* - {a b) {Xhet h)* 



Aa 



{a^){b^) - {aby 

a?){Xhet hb)* - {ab){Xhet ha)* 



(C12) 



(a2)(62) _ i^ab)' 

In the case of the ^-statistic it is useful to introduce a complex amplitude A 

A = Ac + iAg, 



(C13) 



where real amplitudes Ac and As are defined by Eqs. 20 of (Jaranowski & Krolak 2010) and 
a complex filter hg 

hg = {hc + iK)e''''^^\ (CM) 
where real filters he and hg are defined by Eqs. 7 of ( Jaranowski fc Kr61ak|[20l0 ). In complex 



notation the ^-statistic assumes the following simple form (cf Eq. 18 of (Jaranowski & Krolak 

Q 



2010)) 



\{Xhet hg)\'^ 



2D 



where 



D={\hg?). 

The estimator of the complex amplitude A is given by 

(Xhet hg) 



A 



D 



(C15) 
(C16) 

(C17) 
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Grubbs' Test 



The Grubbs' test (Grubbs 1969) is used to detect outliers in a univariate dataset. 

This outlier is removed from the dataset and 



Grubbs' test detects one outlier at a time 
the test is iterated until no outliers are detected 



Grubbs' test is a test of the null hypothesis: 



HO: There are no outliers in the dataset Xi. 
against the alternate hypotheses: 



HI: There is at least one outlier in the dataset Xi. 



Grubbs' test assumes that the data can be reasonably approximated by a normal dis- 
tribution. 

The Grubbs test statistic is the largest absolute deviation from the sample mean in 
units of the sample standard deviation and it is defined as: 

Q _ max|xt — ^1 (1)1) 
a 

where and a denote the sample mean and standard deviation, respectively. 
The hypothesis of no outliers is rejected if 



t2 

a/{2n),n-2 (D2) 



2 + ^a/(2n),n-2 



with ta/{2n),n-2 denoting the critical value of the t-distribution with n — 2 degrees of freedom 
and a significance level of a/{2n). 

We have applied the Grubbs' test to the coarse heterodyne data before analyzing them 
with J-" and Q statistics. We have applied the test to segments of 2^^ data points and we have 
assumed false alarm probability of 0.1%. This resulted in identification of 13 844 outliers 
from the original dataset containing 12 403 138 points. We replaced these outliers with zeros. 
The time series before and after the removal of the outliers are presented in Fig. [6j The 
number of outliers constitutes 0.1% of the total data points in input data resulting in a 
negligible loss of signal-to-noise ratio of any continuous signal present in the data. With 
different methods to identify the outliers used by other searches the number of outliers was 
similar. 
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Fig. 6. — The coarse heterodyne data before (blue) and after (red) removal of the outhers 
using the Grubbs test. The top panel presents the real part of the data and the bottom 
panel the imaginary one. Not all the input data are shown because some outliers are large. 
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of the data respectively. 
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